Absence of Detailed Balance in Ecology 
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Living systems are typically characterized by irreversible processes. A condition equivalent to the reversibility 
is the detailed balance, whose absence is an obstacle for analytically solving ecological models. We revisit 
a promising model with an elegant field-theoretic analytic solution and show that the theoretical analysis is 
invalid because of an implicit assumption of detailed balance. A signature of the difficulties is evident in the 
inconsistencies appearing in the many-point correlation functions and in the analytical formula for the species 
area relationship. 



Spatially explicit models of ecology amenable to analytic solution are rare HKH. This is because ecological measures such 
as the species area relationship (SAR), the relation between the average number of observed species and the sampled area, 
crucially depend on the behavior of the correlation functions at all orders, and any truncation inevitably impairs the results. As a 
consequence, one needs to solve the model in full generality, a task that is highly non-trivial because stochastic theories defined 
on space often have stationary states for which detailed balance (DB) does not hold. 

The main difference between equilibrium stationary states (those for which the DB condition holds) and non-equilibrium sta- 
tionary states is reversibility. The pathways among configurations in the stationary state of an equilibrium system are symmetric 
under time reversal symmetry. In the stationary state of a non-equilibrium system, instead, it is possible to identify a direction 
of time: the 'macroscopic' history of the system with time running backwards or forwards is completely different. For instance, 
if DB is satisfied, the spatio-temporal patterns of declining as well as increasing species' populations should be (on average) the 
same, because under such a condition the decline of a species' population with time running forward is equivalent to an increase 
of the same population, but with the time running backward. However, there are observational findings |4| which show that 
declining species typically have sparse distributions, whereas growing ones are more aggregated. Such patterns emerge because 
the underlying communities are not in equilibrium and call for non-equilibrium models. Such frameworks are usually more real- 
istic but much harder to study. Recently, O'Dwyer and Green ||5l proposed a spatially explicit model (the OG model) for which 
the DB condition is not satisfied. Thus this model is able to show interesting irreversible phenomena, potentially in agreement 
with observed data. For such a model, a very interesting pattern to study is the SAR, for which there is no analytic solution 
based on a microscopic model incorporating fundamental processes. Remarkably, OG were able to derive an explicit formula 
for the SAR by using field theoretical techniques. However, their calculations were based on the incoiTect assumption that DB 
was satisfied. The assumption is not merely a useful (and otherwise harmless) approximation, because such a simplification can 
result in meaningless quantities such as negative probabilities or asymmetric correlation functions. 

A generic Markov process (such as the OG model) is defined via its master equation which specifies the time evolution of 
the probability of any configuration of the system in terms of the transition rates between any two pairs of its microscopic 
configurations: 



where C represents a configuration of the system (e.g., in the OG model this is the set of the numbers of individuals that are 
present in every lattice site) and P{C; t) the probability to find the system in the configuration C at time t. The functions 
W[C C] are the transition rates from the configuration C to the configuration C" (e.g., the buth, death and dispersal rate in 
the OG model). 

An important property of a system is the stationary distribution of its configurations: this is given by the probability, Ps{C), 
which makes the right hand side of eq. [T| equal to zero. Thus, finding Ps{C) is equivalent to calculating the solution of the 
following equation R{C, C) = 0. For some simple systems, the solution can be obtained by assuming that each term in 
the sum is zero, i.e. i?(C, C") = 0. This equation defines the well-known DB condition. However, in general, the stationary 
distribution of a system is described by a probability function that does not satisfy the equation W[G' — )■ C]Ps{C') = W[C ^• 
C']Ps{C) for every C and C . A process whose stationary state obeys DB is said to be an equilibrium stationary state. In this 
case, the stationary state does not exibit any stationary current or flow (e.g. of energy or of particles) and it is possible to prove 
that the stationary solution obey thermodynamic principles and can be described by the poweiful tools of equilibrium statistical 
mechanics |6|. 

Unfortunately, we are not usually given the stationary distribution Ps{C), but only the transition rates W[C C] are known. 
Thus, we need to find out Ps{C) from the transition rates without assuming that the DB condition holds a priori; this is usually 




(1) 



R(C,C') 



2 



a very challenging task. However, it is easy to prove 16] |7] the following equivalent condition known as Kolmogorov's criterion: 
if the DB condition is satisfied, then for every set of configurations {Ci, C2, . . . , C„}: 

W[Ci~^C2]W[C2'^C^\...W[Cn^Ci] = 

W[Cr. ^ C„_i] . . . W[C2 ^ Ci]M^[Ci ^ C„] . 

This equation is a condition that depends only on the transition rates and it always holds for equilibrium steady states. It is also 
a time reversal symmetry property: given a configuration of the system, the probability to come back to the same configuration 
following a path in the space of configurations does not depend on the orientation of the path itself. Thus, eq.|2]suggests a simple 
scheme for assessing whether the OG model has non-equilibrium characteristics: we only need to find a set of configurations 
which do not satisfy the condition in eq. |2] 

The OG model |5| is a birth-death-dispersal neutral process fSl] defined on a regular lattice A (see fig.[T]). The configurations 
of the model are represented by the number of individuals, which are present on each site. The dynamics of different species are 
totally decoupled, i.e., there are no interactions between individuals of different species Il9l 4in . At any given time t, a species 
has a total number of individuals n{t) and any individual can die at a rate d or survive and produce one offspring at a rate b 
if there is at least one individual in the lattice, i.e., 71 > 0. When the species population is zero everywhere, i.e., ?i = 0, a 
speciation/immigration event can occur with speciation rate per site equal to v. The non-trivial spatial dependence of the model 
is given by the birth event: when an individual in the site i produces an offspring, the new individual remains in the same site 
with probability 7 or, alternatively, is located within one of the other sites with probability 1 — 7. The model can be formulated 
with a general dispersion kernel, that defines the probability that the new offspring is located in a site j given a parent in a site 
i. The model was solved only in the case of dispersion to nearest neighbor sites and we focus our analysis on this case. All the 
consequences we present are also valid in the more general case. Thus, the dynamics of the model is completely specified by the 
following transition rates which define the master equation via eq. [T] 

1 - 

. .. . . . ,n|A|} ^- {ni,. . . ,n, + 1,. . . ,n|A|}] = &7nj + ^ + H "^"j'O 

W[{ni, ...,ni,.. . ,n|A|} {ui, ...,n,-l,.. . ,n|A|}] = dn^, 

where |A| is the total number of sites in the lattice and is the population on the site i. As we said, it is sufficient to find a set 
of configuration which do not satisfy the equation |2] to show that the DB does not hold. Figure |2] shows a set of configuration 
which, given the rates defined in eq.|3] does not respect the reversibility condition of eq.|2]and therefore the DB condition is not 
valid for the model. When 7 = 1 there is not dispersion in space, and the model becomes a simple birth-death process for which 
DB holds. 

If a model satisfies DB, the spatio-temporal evolution following a speciation event with time running forwards can be in- 
terpreted (on average) as the spatio-temporal evolution preceding an extinction event with time running backwards. Thus the 
processes after speciation and before extinction events are temporally symmetric. For the OG model, instead, speciation and 
extinction are temporally asymmetric processes (as is observed in nature): the spatio-temporal behavior of species that follows 
a speciation event is completely different (on average) from the spatio-temporal behavior that precedes an extinction event. 

By introducing the moment generator function Z, defined as 

Z{h,t)^ P(ni,...,n|A|;t)e5:..A".''. , (4) 

"1,---,"|A| 

we obtain the following evolution equation: 



dt 



(5) 



where 9{t) is defined as vP{ni = 0, n2 = 0, . . . , n|A| ~ 0; t)/h and /i is the number of neighbors in the lattice (e. g. in a I? 
dimensional square lattice /i — 2D). We can introduce in the previous expression the discrete Laplacian operator Aa defined 
over a generic regular lattice A with periodic boundary conditions 



(Aa/),= ifi~fd^ (6) 

j:|i-«l = l 



obtaining 

dZ{...,h,,...,t) 

dt 
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where cr^ = (1 — 7)//x. In the continuum Hmit (see the Supplementary Materials of 12]), the discrete Laplacian Aa becomes 
the Laplacian operator A and thus we obtain 



J An L SH{x) 6H{x) 



dt 



(8) 



where Aq is the area (volume) of the entire system. 

One can show that the solution proposed by the authors in ref. is equivalent to imposing DB. This condition for the OG 
model (see fig. [TJ can be written as 



d(n, + l)F,^«(...,n, + l,...,n,,...) 



j:\j-i\ = l j 

for any choice of i. The notation P^^ stands for the stationary probability (i.e. the solution of the right hand side of master 
equation set equal to zero) obtained by imposing the DB condition. By introducing into this expression the definition of Z, we 
obtain: 

^ f^.dZf ^ dZf ^ 2 V- 9Zf 

y.\j^i\ = \ 

where 9 is now equal to vP^^ (n\ = 0, . . . , n| a| =0) jh. This expression, in the continuum limit, becomes 

This is exactly the equation solved in ref. jSl . Thus we have shown that the solution obtained by the authors by considering 
eq. 1 1 instead of the right side of eq [8] corresponds to imposing the detailed balance condition which is not vaUd for the OG 



model (as shown in fig.|2]l. 

Assuming DB to hold when it is violated produces incorrect results. An example arises in the Z-point correlation function. 
OG obtained an equation for the moment generating function Z defined in eq. |4j by expanding it, one can obtain the equation 
for the /-points correlation function, {ni-^rii^ . . .rii^') := Gi{ii, . . . , ii). Because of its definition, the correlation function must 
be symmetric under the exchange of its arguments, e.g., for the 3-point correlation function G3(ii, «2, is) = 03(^2, ii, is). 



We can show that applying DB to the OG model (i.e. by solving equation 111, one obtains asymmetric correlation functions. 
Indeed the 3-point correlation function G3(ii, 12, is) does not possess the required symmetry. To see this in an explicit way, 
a useful procedure to simplify and solve the involved differential equations is to apply the Fourier transform and consider the 
function Gs{p_^,p^,p^). Because of the properties of the Fourier transform, we know that the function is symmetric in 

its arguments if and only if G3 is symmetric. By considering the first three correlation functions, we obtain the following 
expression: 

hf) 



-dGi{p^+p^) -d{n) 

(12) 



G2{P,,P,) -^baY~+ib-d) = -ba^pl + ib-d)^^^^ +22) 



^^iPvlrls) -baY, + ib-d) " 

2d{n) - ba^f + {b-d) 
= {baY,~b + d)^ -^(^i +^2+^3) ■ 

The three point correlation function is not symmetric in its arguments. It clearly shows that, on assuming DB in the OG model, 
we obtain inconsistencies in the solution. This inconsistence is valid for all the /-point correlation functions, except the 2-point 
correlation function. 

We have also simulated the model in a simple case: we show in fig.[3]that the analytical results do not match the numerical 
simulations of the model. We consider the OG with just two sites. It is the simplest case in which DB is not valid (i.e. it is 
possible to construct a cycle of configurations as in fig. |2]i. We find a stationary probability by applying DB and show that this 
condition produces inconsistencies. 
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Consider the generic configuration {n, m). The DB conditions can be expressed as 

('dP^^(l,0) = ^P^^(0,0) 

\d{n + l)P°^{n + 1, to) = 6[7?i + (1 - 7)m] P'°'^(n, m) . 
We can thus simply construct a solution by imposing DB. Specifically, we obtain for to and n bigger than zero 



(13) 
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where r(-) is the Gamma function and x = h/d. This expression is clearly inconsistent because it is not symmetric under the 
exchange of its arguments (n, m). We can see that this expression turns to be consistent for 7 = 1 and 7 = 1/2, exactly the 
choices of parameters for which the reversal symmetry in fig. |2]is restored. This inconsistency is shown in fig. |3] where we 
compare the analytical result with the numerical simulations. The case 7 = 1 (which corresponds to = 0) is trivial because 
in this case we are not considering diffusion in space and the sites are independent (see fig. [T]) |9|, and so the DB condition is 
valid. In fact the 3-point correlation function in eq. 12 is symmetric in the case 7=1. When we consider 7 = 1/2 we expect 



the DB not to be valid, because the 3-point correlation function is still non-symmetric. The fact that the DB seems to be valid in 
the simulation of fig.[3]is simply due to the fact that we are considering a very simple case with only two sites. 

We have shown that the solution obtained by applying the DB is not the real stationary solution of the model and that one 
obtains incorrect results. It is not possible to treat that solution as an approximation to the correct one, because there is no 
parameter which is able to quantify the goodness of the approximation. In addition, the 3-point correlation function obtained on 
assuming DB is not symmetric and this leads to incorrect results, regardless of the nature of the approximation. 

One may ask whether the formula obtained by applying the DB could provide a good fit to the data. Unfortunately this is not 
the case. It is known that, at very small areas, the S AR grows linearly with the area ISl [T2 and, in particular, it is equal to the 
number of sampled individuals. By expanding at small areas (r ^ a) the solution obtained by imposing the DB, we obtain: 



S{r) ~ N{r) log( 



ar^ (l + a)(log(5)-7i;) 



(15) 



where 7^ is the Euler-Mascheroni constant (75 w 0.577 . . . ) and N{r) 



is the number of individuals sampled. Note that 



this expression is not close to the expected linear growth S{r) ^ pixr ~ N{r) (where p is the density of individuals), because 
the logarithm is much larger than one (cr/r ^ 1). Moreover, the expansion leads to the meaningless result that the number of 
species sampled is larger than the number of individuals. This wrong prediction is one of the consequences of the use of the DB 
condition. 

The OG model was recently proposed ||5] to describe the stationary properties of an ecosystem. This model appeared to be 
the first solvable model able to predict the empirical Species-Area Relationship, which is one of the most important stationary 
quantities in ecology. OG solved the model in an elegant manner via the moment generator functional. We have shown that 
the OG model is not an equilibrium model, as demonstrated in fig. [2] The OG model is a birth-death-migration process: 
a species changes its population and the individuals diffuse in space, a process which is not reversible. Furthermore, there are 
inconsistencies in the solution on wrongly assuming detailed balance such as the correlation functions no longer being symmetric 
in their arguments. The basic message is that the behaviors and the properties of the stationary state of equilibrium and non- 
equilibrium systems can be distinct. Life is not an equilibrium system nor it is reversible. Developing techniques for studying 
the behavior of models, which do not obey detailed balance, is a necessary but daunting step to understand the physics of non- 
equilibrium processes. An alternative to the OG model based on the Poisson cluster processes, which is solvable and leads to an 
analytic expression for the SAR, has been recently developed 1121 . 
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FIG. 1: Definition of the OG model. The model is non-interacting (different species are totally decoupled) and neutral (the 
parameters of the model do not depend on the species we are considering). Because of these assumptions it is possible to 
consider one species at a time. In the figure, a one dimensional lattice is shown for simplicity. A circle corresponds to a single 
individual and its color represents the species it belongs to. There is no limit to the number of individuals that can be placed in 
each site. The dynamics is defined in terms of three possible main moves: birth, death and speciation. The per capita death rate 

is equal to d, whereas the speciation rate (defined per unit of lattice spacing) is v. The total birth rate is equal to h\ with a 
probability 7, the offspring is placed in the same site as the parent, while, with a probability 1 — 7, it is placed in a different site 
chosen uniformly among the neighboring sites (which corresponds to a dispersal move). 
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b^d^[Yn+(l-Y)m][Ym+(l-Y)(n + l)][m+l][n + l] 
b2d2[Ym+(l-Y)n][Yn+(l-Y)(m+l)][m + l][n + l] 

different! 

if Y^Vzj y^l or mxn 

FIG. 2: Detailed Balance in the OG model. In general, detailed balance is not satisfied if the probability to go along a closed 
path in the space of configurations depends on the direction that one chooses (see eq.|2]l. Here we present the simplest 
counterexample to eq.|2]for the OG model. The probability to find the system in any of the four configurations in red and linked 
by the cycle depends on whether one follows the green or blue path. When 7 = 1, detailed balance is valid, whereas, for 
J — 1/2, it is generally not valid (i.e. even though the condition of eq.[2]is valid for this particular cycle, it is not valid 

generally). 
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FIG. 3: Comparison between the simulations of the OG model and the analytical solution obtained by applying the DB 
condition. The blue circles are simulations of the cumulative stationary probabilities Pi {n) = X^feLn X]m=o ^(^i '^)- The red 

triangles represents P^^^{n), obtained by summing over m eq.[l4] which is calculated by imposing DB. As expected, the 
solutions are different, except in the case 7 = 1/2 (as expected from fig.|2]). For this value of 7 the DB is still not valid (3-point 
correlation function is still non-symmetric). This special behavior of the model for 7 = 0.5 is therefore caused by the smallness 
of the system considered in this case, and not because DB is restored. The model is simulated with the following parameters: 

b = 0.9, rf = 1. and = 0.1. 



